#include "ecgdetection.h"
#include "ecgfilter.h"
#include <algorithm>
#include <numeric>

#define LOWER_HZ (10.0)
#define UPPER_HZ (30.0)
#define SUM_FILTER_DELAY_SEC (0.06)
#define WINDOW_SEC (0.090)//(0.160)
#define MIN_RR_SEC (0.220)//(0.3)//(0.4)//(0.200)
#define MAX_RR_SEC (2.0)

EcgDetection::EcgDetection(QObject *parent)
    : QObject{parent}
{}

void EcgDetection::recognition(QList<int>& ecgSignalData)
{
    //低通高通滤波
    EcgFilter ecgFilter;
    for (int i = 0; i < ecgSignalData.size(); ++i) {
        int fdatum = ecgFilter.lp(ecgSignalData[i]);
        fdatum = ecgFilter.hp(fdatum);
        ecgSignalData[i] = fdatum;
    }
    //数据归一化
    QList<float> norSignalData;
    norSignalData.fill(0,ecgSignalData.size());
    normalization(ecgSignalData,norSignalData);
    int a = 0;
    QRS_message qrs_message;
    QRSCalibration(norSignalData, &qrs_message);
    QList<int> qrs;
    for (int i = 0; i < result.size(); ++i) {
        if(result[i] == 1){
            qrs.append(i);
            qDebug() << i;
        }
    }
    qDebug() << qrs;

}

void EcgDetection::normalization(const QList<int> &ecgSignalData,QList<float>& norSignalData)
{

    float maxValue, minValue,avgAvlue;
    // 获取最大值
    // maxValue = *std::max_element(ecgSignalData.begin(), ecgSignalData.end());

    // // 获取最小值
    // minValue = *std::min_element(ecgSignalData.begin(), ecgSignalData.end());

    // 获取平均值
    double sum = std::accumulate(ecgSignalData.begin(), ecgSignalData.end(), 0);
    avgAvlue = sum / ecgSignalData.size();
    for (int i = 0; i < ecgSignalData.size(); i++)
    {
        norSignalData[i] = ecgSignalData[i]  - avgAvlue;
    }
}

int EcgDetection::QRSCalibration(const QList<float> &norSignalData,QRS_message* qrs_message)
{
    int sampleRate = 125;
    int length = norSignalData.size();
    result.fill(0,norSignalData.size());
    int16_t rpos1[20];
    DetectQrsPeaks(norSignalData, result, sampleRate);
    int16_t n = 0;
    for (int16_t i = 0; i < length; i++)
    {
        if (result[i] == 1 && i > 125)//300/4
        {
            if (i - n < 25)
            {
                result[i] = 0;
            }
            n = i;
        }

    }
    return 0;
    //误判检查
    int16_t rpos_count1t = 0;

    for (int16_t i = 0; i < length; i++)
    {
        if (result[i] == 1 && i > 125)//300/4
        {
            if (rpos_count1t > 19)
            {
                ////break;
                return 0;
            }
            rpos1[rpos_count1t] = i;
            qrs_message->rpos2[rpos_count1t] = i;
            rpos_count1t++;
        }

    }
    for (int16_t i = 0; i < rpos_count1t - 3; i++)
    {
        if (abs(rpos1[i + 2] * 2 - rpos1[i] - rpos1[i + 3]) < 10 || abs(rpos1[i + 2] * 2 - rpos1[i] - rpos1[i + 3]) == 17)
        {
            rpos1[i + 1] = 0;
        }

    }
    //误检整理
    int16_t rpos_countt = 0;
    for (int16_t i = 0; i < rpos_count1t; i++)
    {
        if (rpos1[i] != 0)//rpos1
        {
            qrs_message->rpos[rpos_countt] = rpos1[i];//rpos1
            rpos_countt++;
        }

    }
    qrs_message->rpos_count = rpos_countt;
    qrs_message->rpos_count1 = rpos_count1t;
    return 0;
}

int EcgDetection::DetectQrsPeaks(const QList<float> &norSignalData,QList<unsigned char>& result,float rate)
{
    const int WINDOW_SIZE = (int)(WINDOW_SEC * rate);
    const int MIN_RR = (int)(MIN_RR_SEC * rate);
    const int MAX_RR = (int)(MAX_RR_SEC * rate);
    int count, delay;
    int size = norSignalData.size();

    QList<float> filSignalData;
    filSignalData.fill(0,size);
    delay = FilterSignal(norSignalData, filSignalData, rate);
    Normalize(filSignalData);

    QList<float> derivative;
    derivative.fill(0,filSignalData.size());
    ComputeDerivative(filSignalData, derivative);


    Normalize(derivative);
    ArrayPow2(derivative);


    WindowIntegration(derivative, filSignalData, WINDOW_SIZE);
    // QString str1 = "";
    // for (int i = 0; i < filSignalData.size(); ++i) {
    //     str1 +="," +QString::number(filSignalData[i]);
    // }
    // qDebug() << str1;


    delay += WINDOW_SIZE / 2;

    count = Thresholding2(filSignalData, MIN_RR, MAX_RR, result);
    SubtractDelay(result, delay);
    return count;
}

int EcgDetection::FilterSignal(const QList<float> &norSignalData,QList<float>& filSignalData, float rate)
{
    int size = norSignalData.size();
    Filter filter = { 0,0,0,0,0,0,0,0,0,0 };
    for (int i = 0; i < size; i++)
    {
        filSignalData[i] = norSignalData[i];
    }

    filter = InitFilter(UPPER_HZ, rate, FT_LOW_PASS);
    FilterData(&filter, filSignalData);


    filter = InitFilter(LOWER_HZ, rate, FT_HIGH_PASS);
    FilterData(&filter, filSignalData);
    return (int)round(SUM_FILTER_DELAY_SEC * rate);
}

Filter EcgDetection::InitFilter(float cutoff_hz, float rate, FilterType type)
{
    float cutoff = cutoff_hz / rate;
    Filter filter;
    ResetFilter(&filter);


    fill(&filter, cutoff, (type == FT_HIGH_PASS));

    return filter;
}

void EcgDetection::ResetFilter(Filter *filter)
{
    filter->x0 = 0;
    filter->x1 = 0;
    filter->x2 = 0;
    filter->y1 = 0;
    filter->y2 = 0;
}

void EcgDetection::fill(Filter *filter, float cutoff, int is_hpf)
{
    const float B = tan(cutoff * 3.141592);
    const float BB = B * B;
    const float S = 1.0 + 1.414213 * B + BB;

    if (is_hpf) {
        filter->b0 = 1.0 / S;
        filter->b1 = -2.0 * filter->b0;
    }
    else {
        filter->b0 = BB / S;
        filter->b1 = 2.0 * filter->b0;
    }
    filter->b2 = filter->b0;
    filter->a1 = 2.0 * (BB - 1.0) / S;
    filter->a2 = (1.0 - 1.414213 * B + BB) / S;
}

void EcgDetection::FilterData(Filter *filter, QList<float> &data)
{
    int i;
    int size = data.size();
    for (i = 0; i < size; ++i) {
        filter->x2 = filter->x1;
        filter->x1 = filter->x0;
        filter->x0 = data[i];
        data[i] = filter->b0 * filter->x0 + filter->b1 * filter->x1 +
                  filter->b2 * filter->x2 - filter->a1 * filter->y1 -
                  filter->a2 * filter->y2;
        filter->y2 = filter->y1;
        filter->y1 = data[i];
    }
}
//排序后取出百分比所在的值
float EcgDetection::percentiles(const QList<float> &data, int percent)
{
    QList<float> sortedList = data;
    // 对列表进行升序排序
    std::sort(sortedList.begin(), sortedList.end());
    return sortedList[(int)((float)sortedList.size() * (float)percent / 100 + 0.5)];
}

void EcgDetection::Normalize(QList<float> &values)
{
    int i;
    int size = values.size();
    float max_value = values[0];

    for (i = 1; i < size; ++i) {
        if (values[i] > max_value) {
            max_value = values[i];
        }
    }

    for (i = 0; i < size; ++i) {
        values[i] /= max_value;
    }
}
void EcgDetection::ComputeDerivative(const QList<float> &filSignalData, QList<float> &derivative)
{
    int i;
    int size = filSignalData.size();
    for (i = 2; i < size - 2; ++i) {
        float value = (-filSignalData[i - 2] - 2 * filSignalData[i - 1] +
                       2 * filSignalData[i + 1] + filSignalData[i + 2]);
        derivative[i] = value / 8.0;
    }
    derivative[0] = derivative[1] = derivative[2];
    derivative[size - 1] = derivative[size - 2] = derivative[size - 3];
}

void EcgDetection::ArrayPow2(QList<float> &signal)
{
    int i;
    int size = signal.size();
    for (i = 0; i < size; ++i) {
        signal[i] *= signal[i];
    }
}

void EcgDetection::WindowIntegration(const QList<float> &signal, QList<float> &output, int window_size)
{
    int i;
    int size = signal.size();
    float value = 0.0;
    //float inv_window_size = 1.0 / window_size;

    for (i = 0; i < size; ++i) {
        int first = i - (window_size - 1);
        //value += signal[index] * inv_window_size;
        value += signal[i] / window_size;
        if (first > 0) {
            //value -= signal[first - 1] * inv_window_size;
            value -= signal[first - 1] / window_size;
        }
        output[i] = value;
    }
}

int EcgDetection::Thresholding2(const QList<float>& integrated, int min_rr_width, int max_rr_width, QList<unsigned char>& result)
{
    int size = integrated.size();
    int i, count, previous, searchback, searchback_end;
    float spki, npki, threshold1, threshold2;

    spki = npki = 0.0;
    count = 0;
    threshold1 = 0.0;
    threshold2 = 0.0;
    previous = 0;
    searchback_end = 0;
    searchback = 0;

    QList<float> peaki_sample;
    float peaki_50 = 0.0;
    int num = 0;
    for (i = 7; i < size - 6; ++i)
    {
        float peaki;
        peaki = integrated[i];

        if (peaki > integrated[i + 1] && peaki > integrated[i - 1] && peaki > integrated[i - 2] && peaki > integrated[i + 2] && peaki > integrated[i - 6] * 1.5  && peaki > integrated[i + 6] * 1.5)//1.5
        {
            peaki_sample.append(peaki);
            num++;
        }

    }
    peaki_50 = percentiles(peaki_sample, 73);//53


    for (i = 4; i < size - 4; ++i) {
        int is_qrs;
        float peaki;
        if (i - previous > max_rr_width && i - searchback_end > max_rr_width) {
            searchback = 1;
            searchback_end = i;
            i = previous + 1;
            continue;
        }
        if (searchback && i == searchback_end) {
            searchback = 0;
            continue;
        }
        peaki = integrated[i];
        if (peaki < integrated[i - 2] || peaki < integrated[i - 4] || peaki <= integrated[i + 2] || peaki <= integrated[i + 4] || peaki < peaki_50 * 0.32)
        {
            continue;
        }
        is_qrs = 0;
        if (searchback) {
            if (peaki > threshold2) {
                spki = 0.750 * spki + 0.250 * peaki;
                is_qrs = 1;
            }
        }
        else if (peaki > threshold1) {
            spki = 0.875 * spki + 0.125 * peaki;
            is_qrs = 1;
        }
        if (is_qrs) {
            if (count == 0 || i - previous >= min_rr_width) {
                result[i] = MARK_QRS;
                ++count;
            }
            else if (integrated[previous] < peaki) {
                result[previous] = MARK_NO_QRS;
                result[i] = MARK_QRS;
            }
            previous = i;
        }
        else {
            npki = 0.875 * npki + 0.125 * peaki;
        }
        threshold1 = npki + 0.25 * (spki - npki);
        threshold2 = 0.5 * threshold1;
        ++i;
    }
    return count;
}

void EcgDetection::SubtractDelay(QList<unsigned char> &qrs_detection_result, int samples_delay)
{
    int i;
    int size = qrs_detection_result.size();

    for (i = samples_delay; i < size; ++i) {
        if (qrs_detection_result[i] == MARK_NO_QRS) {
            continue;
        }
        qrs_detection_result[i] = MARK_NO_QRS;
        qrs_detection_result[i - samples_delay] = MARK_QRS;
    }
}

